f3 and life-history traits

Row

Body Size

Trophic level

Propagule size

Fecundity

Age at first maturity

Lifespan

Adult lifespan

Pelagic larval duration

Model selection

# --------------------------------------------------------------------------------------------------
# Initial Model:

Call:
glm(formula = as.formula(paste(Y, paste(in.variable, collapse = "+"), 
    sep = "~")), family = binomial(logit), data = data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.222  -1.222   1.133   1.133   1.133  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.1054     0.4595   0.229    0.819

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 26.287  on 18  degrees of freedom
Residual deviance: 26.287  on 18  degrees of freedom
AIC: 28.287

Number of Fisher Scoring iterations: 3

# -------------------------------------------------------------------------------------------------- 
### iter num = 1, Forward Selection by LR Test: + Adult_Lifespan 

Call:
glm(formula = admixture_f3_negative_number ~ Adult_Lifespan, 
    family = binomial(logit), data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6161  -0.6088   0.0166   0.3955   1.6733  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)  
(Intercept)     -5.3303     2.5811  -2.065   0.0389 *
Adult_Lifespan   1.0534     0.5266   2.000   0.0455 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 26.287  on 18  degrees of freedom
Residual deviance: 14.982  on 17  degrees of freedom
AIC: 18.982

Number of Fisher Scoring iterations: 6

--------------- Variance Inflating Factor (VIF) --------------- 
Multicollinearity Problem: Variance Inflating Factor (VIF) is bigger than 10 (Continuous Variable) or is bigger than 2.5 (Categorical Variable)
# ================================================================================================== 
*** Stepwise Final Model (in.lr.test: sle = 0.15; out.lr.test: sls = 0.15): 

Call:
glm(formula = admixture_f3_negative_number ~ Adult_Lifespan, 
    family = binomial(logit), data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6161  -0.6088   0.0166   0.3955   1.6733  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)  
(Intercept)     -5.3303     2.5811  -2.065   0.0389 *
Adult_Lifespan   1.0534     0.5266   2.000   0.0455 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 26.287  on 18  degrees of freedom
Residual deviance: 14.982  on 17  degrees of freedom
AIC: 18.982

Number of Fisher Scoring iterations: 6

--------------- Variance Inflating Factor (VIF) --------------- 
Multicollinearity Problem: Variance Inflating Factor (VIF) is bigger than 10 (Continuous Variable) or is bigger than 2.5 (Categorical Variable)

Packages used

      My.stepwise           ggrepel            ggpubr          reshape2           stringr           viridis 
          "0.1.0"           "0.9.1"           "0.4.0"           "1.4.4"           "1.4.1"           "0.6.2" 
      viridisLite           R.utils              R.oo       R.methodsS3        kableExtra       formattable 
          "0.4.1"          "2.11.0"          "1.24.0"           "1.8.1"           "1.3.4"           "0.2.1" 
              png            readxl            plotly             dplyr                DT           leaflet 
          "0.1-7"           "1.3.1"          "4.10.0"          "1.0.10"            "0.20"           "2.1.1" 
             ggsn           ggplot2 rnaturalearthdata     rnaturalearth                sf      RColorBrewer 
          "0.5.0"           "3.4.0"           "0.1.0"           "0.1.0"           "1.0-2"           "1.1-3" 
        rmarkdown     flexdashboard 
           "2.11"           "0.6.0" 
---
title: "f3 and life-history traits"
author: "Pierre Barry"
date: "`r format(Sys.time(), '%d %B, %Y, %H:%M')`"
output: 
  flexdashboard::flex_dashboard:
    theme: paper
    orientation: rows
    social: menu
    source_code: embed
    vertical_layout: scroll
---

```{r global, include=FALSE}
list.of.packages <- c("readxl","ggplot2","ggpubr","ggrepel","plotly","My.stepwise")

for (i in list.of.packages){
  if (i %in% installed.packages()[,"Package"] == FALSE){
    install.packages(i)
  }
  eval(bquote(library(.(i))))
}
```

f3 and life-history traits {data-icon="ion-ios-pulse"}
=======================================================================

```{r}
data = read.table(paste(path,"/output/pop_genomics/f3.csv",sep=""),header = TRUE, sep=",",dec=".")
data = data[,-1]

lfh = as.data.frame(read_excel(paste(path,"/data/lfh/lfh.xlsx",sep="")),sheet=1)
lfh = lfh[,-1]
colnames(lfh)[2] = "Species"

data = merge(lfh,data,on=c("Species"))

fst = read.table(paste(path,"/output/pop_genomics/fst.csv",sep=""),header = T, sep=",")
fst = fst[,-1]
fst = fst[order(fst$Species), ]

data = merge(data,fst,on=c("Species"))

data$italic_species = italic_species=c("italic('A. boyeri')",
                                       "italic('A. fallax')",
                                       "italic('C. galerita')",
                                       "italic('C. julis')",
                                       "italic('D. labrax')",
                                       "italic('D. puntazzo')",
                                       "italic('G. niger')",
                                       "italic('H. guttulatus')",
                                       "italic('L. budegassa')",
                                       "italic('L. mormyrus')",
                                       "italic('M. merluccius')",
                                       "italic('M. surmuletus')",
                                       "italic('P. erythrinus')",
                                       "italic('S. cabrilla')",
                                       "italic('S. cantharus')",
                                       "italic('S. cinereus')",
                                       "italic('S. pilchardus')",
                                       "italic('S. sarda')",
                                       "italic('S. typhle')")

admixture_f3_significant = c()
admixture_f3_negative = c()
admixture_f3_negative_number = c()

for (i in 1:nrow(data)){
  for (j in 1:ncol(data)){
    if (is.na(data[i,j])==TRUE){
      data[i,j] = 0
    }
  }
}

for (i in 1:nrow(data)){
  
  if (data$f3_MU_LI_GA[i] < 0 | data$f3_FA_LI_GA[i] < 0){
    admixture_f3_negative[i]="Yes"
    admixture_f3_negative_number[i] = 1
  } else {
    admixture_f3_negative[i]="No"
    admixture_f3_negative_number[i] = 0
  }
  
  if (data$f3_MU_LI_GA_Z[i] < -3 | data$f3_FA_LI_GA_Z[i] < -3){
    admixture_f3_significant[i]="Yes"
  } else {
    admixture_f3_significant[i]="No"
  }
  
}

data$admixture_f3_negative = admixture_f3_negative
data$admixture_f3_negative_number = admixture_f3_negative_number
data$admixture_f3_significant = admixture_f3_significant

```

Row {.tabset data-height=1000}
-----------------------------------------------------------------------

```{r}
plot_legend =c("Body size (cm)",
               "Trophic level",
               "log[Fecundity (eggs/day)]",
               "Propagule size (mm)",
               "Age at first maturity (years)",
               "Lifespan (years)",
               "Adult lifespan (years)",
               "Pelagic Larval Duration (days)",
               "Hudson's Fst - Gulf of Lion - Costa Calida",
               "Hudson's Fst - Gulf of Lion - Algarve",
               "Hudson's Fst - Gulf of Lion - Bay of Biscay",
               "Hudson's Fst - Costa Calida - Algarve",
               "Hudson's Fst - Costa Calida - Bay of Biscay",
               "Hudson's Fst - Algarve - Bay of Biscay"
               )

my_comparisons <- list( c("No","Yes"))
```

### Body Size

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Body_Size,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Body size (cm) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Body size (cm)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()

ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Trophic level

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Trophic_Level,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Trophic level : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Trophic level") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Propagule size

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Propagule_Size,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Propagule size (cm) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Propagule size (mm)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Fecundity

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Fec,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]
data_tmp$y = log(as.numeric(data_tmp$y))

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Fecundity (eggs/day): ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("log[Fecundity (eggs/day)]") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Age at first maturity

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Age_Mat,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Age at first maturity (years) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Age at first maturity (years)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Lifespan

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Lifespan,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Lifespan (years): ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Lifespan (years)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Adult lifespan

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$Adult_Lifespan,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Adult lifespan (years) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Adult lifespan (years)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Pelagic larval duration

```{r}
data_tmp = data.frame(x=data$admixture_f3_negative,
                      x1=data$admixture_f3_negative_number,
                        y=data$PLD,
                        Species = data$Species_plot,
                        italic_species = italic_species)
  
data_tmp = data_tmp[data_tmp$y!="NA",]

model<-glm(x1 ~ y,
          data = data_tmp,
          family=binomial)

p<-ggplot(data_tmp,aes(x = x, y = y)) +
  geom_violin(alpha = 0.1,aes(fill=x),col="white") +
  geom_jitter(aes(text=paste("Admixture: ", x, "\n Pelagic larval duration (days) : ",y,"\n Species : ",data_tmp$Species,sep="")),size = 3, alpha = 0.5, width = 0.1,col="black") +
  theme_classic() +
  #labs_pubr() +
  xlab("Presence of admixture indicated by negative f3") +
  ylab("Pelagic Larval Duration (days)") +
  geom_text_repel(label=data_tmp$italic_species,col="black",parse=T) + 
  stat_compare_means(comparisons = my_comparisons,
                     method = "t.test",
                     label = "p.format") +
  scale_fill_viridis_d(begin = 0.25,end=0.75) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  ggtitle(paste("p = ",round(coefficients(summary(model))[2,4],4),sep="")) +
  coord_flip()
  
ppp<-ggplotly(p,height=500,tooltip=c("text"))  %>% partial_bundle() 
ppp
``` 

### Model selection

```{r}
My.stepwise.glm(Y = "admixture_f3_negative_number", 
                variable.list = c("Body_Size","Trophic_Level","Propagule_Size",
                                  "Age_Mat","Lifespan","Adult_Lifespan",
                                  "Hermaphrodism","Parental_Care","PLD"
                                  #,"FST_HUDSON_LI_MU","FST_HUDSON_LI_FA","FST_HUDSON_LI_GA",
                                  #"FST_HUDSON_MU_FA","FST_HUDSON_MU_GA","FST_HUDSON_FA_GA"
                                  ),
                data = data, myfamily = "binomial")
``` 

                
Packages used {data-icon="fa-map"}
=======================================================================

```{r}
installed.packages()[names(sessionInfo()$otherPkgs), "Version"]
```